# Analyzing metabolomics data for Mutographs ----
# 03_sensitivity_analysis
# Testing sensitivity of models by shuffling the rows to create "random" metabolites

rm(list = ls())

# 0 - Definition of libraries, paths and functions ----

Sys.setlocale("LC_TIME", "C")

library(stringr)
library(tidyr)
library(Hmisc)
library(MASS)
library(data.table)
library(dplyr)
library(forcats)
library(tibble)
library(openxlsx)
library(ggplot2)
library(cowplot)
library(scales)
library(patchwork)
library(ggpubr)

# Loading previously processed data
metabo_RAW <- read.csv("data/metabolomics_normalized_data.csv", check.names=FALSE)
metabo_PROC <- readRDS("output/metabo_PROC")
metabo_MAN <- readRDS("output/metabo_MAN")
metabo_pearson <- readRDS("output/metabo_pearson")
proxy_table <- readRDS("output/proxy_table")

features <- names(metabo_RAW)[str_detect(names(metabo_RAW),"@")]
metabolites <- names(metabo_PROC)[str_detect(names(metabo_PROC),"@")]

cosmic_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_COSMIC_SBS96_abs_mutations.csv") %>% 
  rename(donor_id=X, SBS1536A_cosmic = SBS1536A, SBS1536B_cosmic = SBS1536B, SBS1536F_cosmic = SBS1536F, SBS1536I_cosmic = SBS1536I)
dbs_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_COSMIC_DBS78_abs_mutations.csv") %>% rename(donor_id=X)
id_attrib <- read.csv("data/sigs/pruned_attribution_RCC_Manuscript_COSMIC_ID83_abs_mutations.csv") %>% rename(donor_id=X)

sigs_cosmic <- names(cosmic_attrib)[2:15]
dbs_cosmic <- names(dbs_attrib)[2:6]
id_cosmic <- names(id_attrib)[2:10]

res_reg_cosm <- readRDS("results/regressions/res_reg_cosm")
res_reg_cosm_sro <- readRDS("results/regressions/res_reg_cosm_sro")

# Creating random metabolites with permutation of the real ones - 10 times (9440 metabolites) ----
set.seed(42)
l_rdm <- list()
for (i in 1:10){
  feature_rdm <- metabo_MAN %>% 
    select(all_of(metabolites)) %>% 
    slice(sample(1:n())) # Shuffling rows
  names(feature_rdm) <- str_c("X", (1+(i-1)*length(feature_rdm)):(i*length(feature_rdm)))
  l_rdm[[i]] <- feature_rdm
}
metabo_RDM <- bind_cols(l_rdm[[1]],l_rdm[[2]],l_rdm[[3]],l_rdm[[4]],l_rdm[[5]],l_rdm[[6]],l_rdm[[7]],l_rdm[[8]],l_rdm[[9]],l_rdm[[10]]) %>% 
  bind_cols(select(metabo_MAN, donor_id, all_of(sigs_cosmic), all_of(dbs_cosmic), all_of(id_cosmic), country, sex, age_diag, bmi, batch, acq_order))
rdm_features <- names(metabo_RDM)[str_starts(names(metabo_RDM),"X[0-9]")]

# Running same regressions for random features with signatures - All cases ----

sparseness_sigcosm <- colSums(select(metabo_MAN, all_of(sigs_cosmic))==0)/NROW(metabo_MAN)*100
reg_cosm_rdm <- list()
for (sig in sigs_cosmic){
  ltmp <- list()
  for (rdm in rdm_features){
    
    df_mod <- metabo_RDM %>% select(rdm, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    
    # Model: Quasipoisson if sparseness below 30%, logistic (A/B median) else
    if (sparseness_sigcosm[sig] >= 30){
      lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    } else {
      lreg <- glm(get(sig) ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "quasipoisson")
    }
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = rdm)
    ltmp[[rdm]] <- resreg
  }
  reg_cosm_rdm[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# Adding logistic for SBS1 and SBS1536A/B (cosmic) as sparseness is < 30%. Others are too zero-inflated (>50%)
for (sig in c("SBS1", "SBS1536A_cosmic", "SBS1536B_cosmic")){
  ltmp <- list()
  for (rdm in rdm_features){
    
    df_mod <- metabo_RDM %>% select(rdm, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = rdm)
    ltmp[[rdm]] <- resreg
  }
  reg_cosm_rdm[[str_c(sig, "_logistic")]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=str_c(sig, "_logistic"), pval_adj = p.adjust(pval, "bonferroni")) 
}
# saveRDS(reg_cosm_rdm, "results/regressions/reg_cosm_rdm"))

res_reg_cosm_rdm <- rbindlist(reg_cosm_rdm) %>% 
  group_by(signature) %>% 
  mutate(sig_label = factor(case_when(signature == "SBS1536B_cosmic" ~ "SBS40a", signature == "SBS1536A_cosmic" ~ "SBS40b", signature == "SBS1536F_cosmic" ~ "SBS40c",
                                      signature == "SBS22" ~ "SBS22a", signature == "SBS1536I_cosmic" ~ "SBS22b", TRUE ~ signature),
                            levels = c("SBS1","SBS2","SBS4","SBS5","SBS12","SBS13","SBS18","SBS21","SBS22a","SBS22b","SBS40a","SBS40b","SBS40c","SBS44")))

# saveRDS(res_reg_cosm_rdm, "results/regressions/res_reg_cosm_rdm")

thresholds_cosm <- res_reg_cosm_rdm %>% 
  group_by(signature) %>% 
  summarize(maxpval = max(-log10(pval))) %>% 
  mutate(bonflim = -log10(0.05/944),
         threshold = pmax(bonflim, maxpval)) %>% 
  left_join(distinct(select(res_reg_cosm, signature, sig_label)))
# saveRDS(thresholds_cosm, "results/regressions/thresholds_cosm")

hits_cosm <- res_reg_cosm %>%
  left_join(thresholds_cosm) %>% 
  filter(-log10(pval) >= threshold) %>% 
  mutate(mass_mz = as.numeric(word(metabo,sep="@"))+1.0078,
         rt = as.numeric(word(metabo,2,sep="@")))
# saveRDS(hits_cosm, "results/regressions/hits_cosm")

# Running regressions - Only in Serbia & Romania for SBS22 and SBS1536I (SBS22b) ----
sparseness_sigcosm_ro <- colSums(select(filter(metabo_MAN, country %in% c("Romania", "Serbia")), all_of(sigs_cosmic))==0)/NROW(filter(metabo_MAN, country %in% c("Romania", "Serbia")))*100
# SBS22 and SBS1536I still super sparse in Romania and Serbia (> 40%), so logistic regression is used
reg_cosm_rdm_sro <- list()
for (sig in c("SBS22", "SBS1536I_cosmic")){
  ltmp <- list()
  for (rdm in rdm_features){
    
    df_mod <- metabo_RDM %>% filter(country %in% c("Romania", "Serbia")) %>% select(rdm, sig, country, age_diag, sex, bmi, batch, acq_order)
    names(df_mod) <- c("metabo", names(df_mod[2:8]))
    df_mod <- df_mod %>% mutate(sig_cat = ifelse(get(sig) > median(get(sig)), 1, 0))
    lreg <- glm(sig_cat ~ metabo + age_diag + sex + bmi + batch/acq_order, data = df_mod, family = "binomial")
    resreg <- coef(summary(lreg))[,c(1,4)] %>% as.data.frame() %>% rownames_to_column()
    names(resreg) <- c("rowname","Estimate","pval")
    resreg <- resreg %>% 
      mutate(metabo = rdm)
    ltmp[[rdm]] <- resreg
  }
  reg_cosm_rdm_sro[[sig]] <- rbindlist(ltmp) %>% filter(rowname == "metabo") %>% mutate(signature=sig, pval_adj = p.adjust(pval, "bonferroni"))
}
# saveRDS(reg_cosm_rdm_sro, "results/regressions/reg_cosm_rdm_sro"))

res_reg_cosm_rdm_sro <- rbindlist(reg_cosm_rdm_sro) %>% 
  group_by(signature) %>% 
  mutate(sig_label = factor(case_when(signature == "SBS1536B_cosmic" ~ "SBS40a", signature == "SBS1536A_cosmic" ~ "SBS40b", signature == "SBS1536F_cosmic" ~ "SBS40c",
                                      signature == "SBS22" ~ "SBS22a", signature == "SBS1536I_cosmic" ~ "SBS22b", TRUE ~ signature),
                            levels = c("SBS1","SBS2","SBS4","SBS5","SBS12","SBS13","SBS18","SBS21","SBS22a","SBS22b","SBS40a","SBS40b","SBS40c","SBS44")))
# saveRDS(res_reg_cosm_rdm_sro, "results/regressions/res_reg_cosm_rdm_sro")

thresholds_cosm_sro <- res_reg_cosm_rdm_sro %>% 
  group_by(signature) %>% 
  summarize(maxpval = max(-log10(pval))) %>% 
  mutate(bonflim = -log10(0.05/944),
         threshold = pmax(bonflim, maxpval)) %>% 
  left_join(distinct(select(res_reg_cosm, signature, sig_label)))
# saveRDS(thresholds_cosm_sro, "results/regressions/thresholds_cosm_sro")

hits_cosm_sro <- res_reg_cosm_sro %>%
  left_join(thresholds_cosm_sro) %>% 
  filter(-log10(pval) >= threshold) %>% 
  mutate(mass_mz = as.numeric(word(metabo,sep="@"))+1.0078,
         rt = as.numeric(word(metabo,2,sep="@")))
# saveRDS(hits_cosm_sro, "results/regressions/hits_cosm_sro")

# Replacing hits all countries by only Romania & Serbia for SBS22a and SBS22b
hits_cosm_all <- hits_cosm %>% 
  filter(!sig_label %in% c("SBS22a","SBS22b")) %>% 
  bind_rows(hits_cosm_sro)
